The criticality of self-assembled rigid rods on triangular lattices 
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The criticality of self-assembled rigid rods on triangular lattices is investigated using Monte Carlo 
simulation. We find a continuous transition between an ordered phase, where the rods are oriented 
along one of the three (equivalent) lattice directions, and a disordered one. We conclude that 
equilibrium polydispersity of the rod lengths does not affect the critical behavior, as we found that 
the criticality is the same as that of monodisperse rods on the same lattice, in contrast with the 
results of recently published work on similar models. 

PACS numbers: 64.60Cn, 61.20.Gy 



Recently experimentalists have acquired the ability to 
control the interactions between colloidal particles, with 
dimensions in the nanometer-to-micrometer range, pro- 
viding new windows into the structural and thermody- 
namic behavior of colloidal suspensions. Of particular 
interest are the so-called "patchy colloids", the surfaces 
of which are patterned so that they attract each other 
via discrete bonding sites (patches) of tunable number, 
size and strength. Some of the collective properties of 
patchy colloids are being intensively studied with the- 
ory and simulations of primitive models and a number of 
results have been obtained 

In systems with two bonding sites per particle, only 
(polydisperse) linear chains form and there is no liquid- 
vapor phase transition |2j]. If the linear chains are suf- 
ficiently stiff, however, they may undergo an ordering 
transition, at fixed concentration, as the temperature de- 
creases below the bonding temperature. The description 
of self-assembled rods has to consider not only the effects 
of equilibrium polydispersity but also the polymerization 
process. In this context, we proposed a model of self- 
assembled rigid rods (SARR), composed of monomers 
with two bonding sites that polymerize reversibly into 
polydisperse chains 0] and carried out extensive Monte 
Carlo (MC) simulations to investigate the nature of the 
ordering transition of the model on the square lattice 
The polydisperse rods undergo a continuous order- 
ing transition that was found to be in the 2D Potts q=2 
(Ising) universality class, as in similar models where the 
rods are monodisperse [5| . This finding refutes the claim 
that equilibrium polydispersity changes the nature of the 
ordering transition of rigid rods on the square lattice from 
Ising to the Potts q=l (percolation) universality class 
and questions a more recent one for a similar model on 



triangular lattices [1, 0| ■ 

In this note we (re)examine the original model of 
SARR [3:, 4] on triangular lattices (TL) using MC simu- 
lations. In the model a site is either empty or occupied 
by one monomer. Each monomer has two interacting 
patches pointing in opposite directions, ±s. The orien- 
tation (state) of the monomers, defined by the direction 
of the bonding patches, is restricted to the (three) lattice 
directions. The interaction potential can be described 
as follows: Provided that two particles, i and j, occupy 
nearest-neighbor (NN) sites and provided that they are in 
the same state, the energy is lowered by an amount e only 
if their orientations arc fully aligned with the lattice vec- 
tor Vij (See Fig. [T]). The bonding energy favors the self- 
assembly of rod- like lattice polymers (straight chains). 

The Grand canonical Hamiltonian of the system is 
given by: 
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i=l k=l 1=1 

(1) 

where i labels a lattice site, dik',k — 1, • • • ,p, are unit 
vectors along the p lattice directions {p — i for TL); s(ri) 
denotes the orientation at a given lattice site: s(ri) = 
for an empty site, while for occupied sites s(ri) is equal 
to one of dfe vectors; M is the total number of sites; 
,f{x) = 1 if a; = 1, and zero otherwise; and /i is the 
chemical potential. 

An ordering transition will occur as the average rod 
length increases. In the ordered phase the rods will align 
preferentially along one of the p lattice directions (See 
Fig. 1). In this model, the only attractions between pairs 
of NN monomers are bonding ones. Additional lateral 
interactions that promote the condensation of monomers. 
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FIG. 1: Examples of ordered (left) and disordered (right) 
configurations for the SARR model on triangular lattices. 
Monomers are represented with thick segments lying on the 
lattice sites. Two nearest-neighbor monomers interact (and 
form a bond) if the corresponding segments are in a head-to- 
tail configuration. 



leading to a competition between ordering of SARR and 
monomer condensation, are not considered, in line with 
the orig inal SARR model Q on the square lattice d, 
0]. The presence of only two bonding patches and the 
absence of lateral interactions strongly suggest that the 
present model does not exhibit a discontinuous liquid- 
vapor transition at low temperatures [H, 0, [![ . 

The investigation of the ordering transition of SARRs 
on the TL is carried out through the analysis of the order 
parameter [Tf, 



(5 = 



(2) 



where Nk is the number of monomers with orientation 
ak- 
in the Grand Canonical Ensemble, at a fixed chemical 
potential, the critical temperature, Tc = Tc{fJ.) is found 
by extrapolating to the thermodynamic limit the finite- 
size pseudo critical temperatures, Tc{L), which may be 
defined in various ways (L being the length of the rhom- 
bic simulation box; M = L^). Given the symmetry of the 
model, one can use the fourth-order cumulant of the or- 
der parameter distribution gi =< 5"^ > / < 5^ >^ , 
to define Tc{L) as the temperature where the finite-size 
cumulant g^iT) takes the universal value, (74cjfor a given 
universality class and boundary conditions 0- We as- 
sume that the criticality of polydisperse rods on the TL 
is the same as that of monodisperse rods on the same 
lattice, i.e.. Potts g = 3 We emphasize that this 
assumption is made for (computational) convenience and 
does not constrain the determination of the critical be- 
havior, as discussed below. 

Although the value of g/^c for the Potts q=2 model 
on the square lattice is well known [isl . we have not 
found in the literature reliable estimates of (74c for the 
Potts (7 = 3 model on TL with periodic boundary con- 
ditions, rhombic boxes, and order parameters defined as 
in Eq.(l2]). We have therefore estimated its value by run- 
ning simulations of the Potts q = 3 model at the critical 
temperaturefisl . with the same box shape and bound- 
ary conditions, using the Swendsen-Wang algorithm [l6j . 



Different system sizes in the range 12 < L < 96 were con- 
sidered. The results were fitted to the scaling equation 



g4{L,Tc) = gic + aL^ 



(3) 



where g^c, «, and yi are obtained from fits of the simu- 
lation results or, alternatively, yi (the critical exponent 
associated to the so-called irrelevant field) is set to the 
theoretical value yi — —4/5 [3, [l^. In the first case we 
find = -0.74 ± 0.10 and g^c = 1-168 ± 0.002; while 
setting yi = —4/5 leads to 54c = 1.167 ± 0.001. We have 
used the latter values in the finite-size scaling analysis 
reported below. 

We carried out coupled Grand Canonical Ensemble MC 
simulations. For a fixed value of fi several values of the 
temperature, T, = To + lAT (with i = 0, 1, • • • , A'^r - 1), 
are sampled in a single MC run using a simulation tem- 
pering algorithm 20] . This is achieved using a probabil- 
ity function given by; 



P(S*^r,) = l](T,)cxp [~H{S^')/kBT, 



(4) 



In order to obtain good sampling over all temperatures 
one has to use an appropriate weight function ^{Ti). 
This was computed through an equilibration procedure 
following the usual strategies of the Wang-Landau-type 
algorithms j20l - [2^ . This simulation temperiiig algorithm 
is known to enhance the sampling efficiency [23|. 

After preliminary runs to locate the critical region we 
run long simulations using typically between Nt = 20 
and Nt — 40 values of Ti around the critical tempera- 
ture. At each fi we considered different system sizes. As 
the interactions are restricted to NN, the lattices are split 
into three sublattices, where the sites do not interact en- 
ergetically. Simulation runs are organized in sweeps. In 
a sweep we update the state of every site and attempt 
one temperature change. This is done by considering se- 
quentially the three sublattices; we select for each site a 
new state (fc = 0, 1, 2, 3) (with k = denoting an empty 
site) with probabilities depending on the interaction en- 
ergy, the value of fj, and of the current temperature. After 
updating all sites we attempt a temperature change by 
choosing at random (with equal probabilities) increas- 
ing or decreasing the current temperature by an amount 
AT, and accept or reject the change by considering the 
probability given by Eq. (H)), and the usual Metropolis 
criterion I9j. The length of a simulation run was 2 x 10® 
sweeps, and the results were split into twenty blocks of 
10^ sweeps for subsequent error analysis. In Figure [2] 
we illustrate the results for the order parameter close to 
the transition temperature, at two values of the chemical 
potential. 

System-size dependent pseudo-critical temperatures, 
TcL = Tc{L) are computed by the matching criterion 
0, 



g4{L,TcL,ti) = gic 



(5) 
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FIG. 2: Results for the order parameter 5 as a function of the 
reduced temperature for different system sizes (as indicated in 
the legends); left panel fi/e — —0.90, right panel p/e = —0.95. 



where for convenience we set g^c to the universal value 
of the Potts q = 3 universality class. Critical tempera- 
tures, Tc, were extrapolated by fitting the values Tc{L) 
to scaling equations of the form, 



Tc{L) =Tc + aL 



-6 



(6) 



In order to avoid biasing the analysis we used two val- 
ues for the exponent b: b = {1 + 9)/v HI with 

— —yi^ and V = 5/6 for Potts g = 3 scaling, 
and b = {l/v)q^i — 0.75 for Potts q = l [15]; v and 9 are 
respectively the correlation length and Wegner's correc- 
tion to scaling exponents. However, the two values of Tc 
were found to be very close. The critical temperatures 
collected in Table H] are those computed using the Potts 
9 = 3 scaling, which are consistent with the temperatures 
where the Binder cumulants cross (see Figure [3]). Notice 
that the crossing of the 54 (T) curves for different values 
of L deviates slightly from the computed value of (74c. In 
order to constrain as little as possible the analysis of the 
criticality of SARRs on the TL we have computed sec- 
ondary error bars (shown between curly brackets in Table 
m that include the estimates for the critical temperature 
found using the Potts q = I scaling. 

The critical behavior of the model was investigated by 
analyzing the system-size dependence of various proper- 
ties at the extrapolated critical temperature. We fit the 
simulation results for a given property at Tc to the ex- 
pected scaling relation [9| 5{L) cx L~^/^ , x{L) oc L'^^", 
and (5 In < S{L) > /dT)^, cx L^/", where x{L) = 
[< (5^ > - < (5 >2] /ksT. (3 and 7 are the critical ex- 
ponents for the order parameter and the susceptibility, 
respectively. In addition, the quantities proportional to 
the second derivatives of the Grand Potential per unit 
volume with respect to the temperature and/or chemical 
potential, {d{H / M) / dT) {dp/dT)^, and {dp/dp)T, are 
fitted to non-linear equations of the form. 



System 


Ai/e = -0.95 


At/e = -0.90 


p oo(p = 1) 


n 


11 


13 


9 




84-192 


72-192 


60-144 


T* 

' c 


0.25336(4){11} 0.29006(4){10} 0.47637(4){19} 


Pc 


0.597(3){9} 


0.688(3){6} 






0.110(22){57} 


0.115(18){43} 


0.126(9){31} 


j/iy 


1.69(7){19} 


1.72(5){12} 


1.70(4){12} 


Ijv 


1.27(8){20} 


1.28(6){12} 


1.21(4){12} 


{alv) 


0.40(28){38} 


0.43(25){33} 


0.45(28){32} 



TABLE I: Finite-size scaling results from simulation: n is the 
number of system sizes used to compute critical properties 
and effective critical exponents. Lmin and Lmai are the mini- 
mum and maximum system sizes used in the finite-size scaling 
analysis. The results shown for ajv were computed from the 
scaling of (9p/9r)^, except for the full lattice case (p = 1) 
where {d[H/M)/dT) was used. For finite p similar values of 
a/u were obtained using {dp/dfi)T or {d{H/M)/dT)fj,. The 
critical exponent ratios for Potts g = 1 universality class are 
P/u = 5/48 ~ 0.104, j/u = 43/24 ~ 1.792, 1/u = 3/4 and 
a/u — —1/2. The corresponding values for Potts q = 3 uni- 
versality class are P/u = 2/15 ~ 0.133, ■y/u = 26/15 ~ 1.733, 
l/u — 6/5, and a/u = 2/5|l^. Error bars are given in paren- 
theses (curly brackets) in units of the last digit of the corre- 
sponding quantity. 




(7) 



FIG. 3: Fourth order Binder cumulant at constant p as a 
function of T, for different system sizes; left panel jj,/e = 
—0.90, right panel p/e = —0.95. Horizontal lines depict the 
estimate of g^c for the two dimensional Potts q=3 universality 
class. Vertical lines delimit the unbiased estimates of the 
critical temperature (See the text for the details). 



where p is the density (fraction of occupied sites), and a is 
the specific heat critical exponent. In Fig. |3]we illustrate 
the p, and of its derivative {dp/dT)i^i at p/e = —0.95 
around the critical temperature. In TableUwe collect the 
estimates for the different critical exponents (or exponent 
ratios) . The uncertainty in the estimate of Tc was taken 
into account and as we did for the critical temperatures, 
two estimates of the error bars are given, with the second 
one corresponding to error bars that are sufficiently large 
to include the critical temperature found using the Potts 
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FIG. 4: Density and derivative of the density with respect to 
the temperature at constant chemical potential /i/e = —0.95. 



q — I scaling. 

For completeness we have computed the critical densi- 
ties, by fitting the results to 24| : 



(8) 



of the effective exponents v, and a/v are not compatible 
with Potts 9=1 critical behavior. The deviations ob- 
served in the crossings of g4{T) for different system sizes 
from the estimated value of g^c is most likely due to the 
importance of scaling corrections (low absolute value of 
Vi)- 

In previously published work ^ we discussed the rea- 
sons for the apparent q = I critical behavior observed 
by Lopez et al. on the square lattice @. The appar- 
ent = 1 behavior observed by the same authors on the 
TL [3| results also from using the density as the scal- 
ing variable. In fact, a simple but revealing analysis by 
Fisher [2^, shows that fixing the density in models such 
as those discussed here, corresponds to introducing a con- 
straint that renormalizes the critical exponents. For the 
Potts (7 = 3 universality class the renormalized correla- 
tion length exponent I'x is t^x = — a) = 5/4, which 
is close to the value of for the q — I universality class 
Vq^i = 4/3, reported by Lopez et al. Q- 
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using the value of z/ for the Potts q = 3 universality class. 

The results in table U clearly indicate that the critical 
behavior of the SARR model on the TL is much bet- 
ter described within the Potts q = 3 than within the 
q = 1 universality class. Given that the former criti- 
cal behavior was observed in the monodisperse case, we 
conclude that equilibrium polydispersity does not affect 
the critical behavior of rigid rod models, in contrast with 
the conclusion of Lopez et al. Q- Even considering the 
largest error bars on the critical temperature, the values 
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